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We review and further develop the recently intro- 
duced numerical approach [Phys. Rev. Lett. 86, 5031, 
(2001)] for scattering calculations based on a so called 
pseudo-time Schrodinger equation, which is in turn a modi- 
fication of the damped Chebyshev polynomial expansion 
scheme [J. Chem. Phys. 103, 2903, (1995)]. The method 
utilizes a special energy-dependent form for the absorb- 
ing potential in the time-independent Schrodinger equa- 
tion, in which the complex energy spectrum is mapped 
inside the unit disk Eu — » Uk, where Uk are the eigen- 
values of some explicitly known sparse matrix U. Most 
importantly for the numerical implementation, all the 
physical eigenvalues ut are the extreme eigenvalues of 
U (i.e., ut ~ 1 for resonances and Uk = 1 for the bound 
states), which allows one to extract these eigenvalues 
very efficiently by harmonic inversion of a pseudo-time au- 
tocorrelation function y(t) — r^U 1 ^ using the filter diago- 
nalization method. The computation of y(t) up to time t — 
2T requires only T sparse real matrix-vector multiplica- 
tions. We describe and compare different schemes, effec- 
tively corresponding to different choices of the energy- 
dependent absorbing potential, and test them numer- 
ically by calculating resonances of the HCO molecule. 
Our numerical tests suggest an optimal scheme that pro- 
vide accurate estimates for most resonance states using 
a single autocorrelation function. 

Key Words: quantum scattering, absorbing potential, 
resonances, pseudo-time Schrodinger equation, Cheby- 
shev polynomial expansion, iterative diagonalization, 
harmonic inversion, filter diagonalization method. 



Introduction. In this paper we present a detailed de- 
scription and further generalization of the methodology de- 



veloped in the preceding works [1, 2, 3, 4, 5] for the efficient 
numerical solution of the quantum scattering problem as- 
sociated with the time-independent Schrodinger equation 



(HiP){r) = Ei/j(r) 



(1) 



Eq. 1 may possess bound states with real energies E 
and wavefunctions ^p(r) exponentially localized in space, 
and resonance states (Siegert states [6]) having complex 
energies with Im E < 0. The latter behave like bound 
states in some compact subset f2 of the configuration space, 
but eventually grow exponentially outside of f2, due to the 
outgoing asymptotic boundary conditions. 

The corresponding bound state problem is conceptually 
simple and very well understood. The numerical solu- 
tion of (1) at dissociation energies is much more difficult 
as it requires the solution of a boundary value problem. 
One can avoid the latter by the use of a so-called optical 
(or absorbing) potential W(r) with Im W(r) < (in the 
sense that i(W — W*) is positive semidefinite, where * de- 
notes conjugate transposition), that vanishes for r G Q and 
smoothly grows outside £1. Numerically, this has negligible 
effect on the scattering solutions ij)(r) of 



(Hip)(r) = (E-W{r))ip(r) 



(2) 



inside f2, i.e., the physically relevant region, and damps 
them outside SI [7]. In other words, the complex absorb- 
ing potential forces the resonance solutions to behave like 
bound states everywhere without significantly affecting the 
energies E. In this framework the physically relevant part 
of the system is, therefore, dissipative and satisfies (1) only 
for r € fl. Moreover, a general multichannel scattering 
problem can be considered with a numerically convenient 
form of W(r), independent of the choice of coordinate sys- 
tem. The price for these benefits is that the originally 
hermitian problem becomes non-hermitian; but it is gen- 
erally still complex symmetric. 
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Although to satisfy Im E < one only needs Im W < 0, 
traditionally one simply uses a negative imaginary poten- 
tial W = —iT and gets the nonhermitian eigenvalue prob- 
lem (H — iT)tp = Erp. The latter is generally much easier 
to handle numerically than the boundary value problem 
(1). As we already saw in [1, 5], energy-dependent choices 
W = We {r) are particularly useful. 

The introduction of the absorbing potential leads to the 
damped Green's function [8, 9, 10] 

G W {E) := (H-E + We)- 1 - (3) 

We assume that under suitable conditions on We{t), sim- 
ilarly to the traditional case W = —iT [10], Gw(E) con- 
verges for any real E (and also for Im E > 0) weakly to 
the ordinary Green's function 

G{E) = \\m{H - E - is)' 1 . 

ej.0 

Practically, one usually needs to evaluate only certain ma- 
trix elements (jr r G(E)ip, the basic numerical objects of 
quantum physics, from which most other quantities of in- 
terest (scattering amplitudes, reaction rates, etc.) can 
be computed (see, e.g., refs. [8, 9, 11, 12]). If both 
and ?p have support in f2, they are well approximated by 

4> t g w (e)tP. 

Unfortunately, for very large systems with high den- 
sity of states one may encounter numerical difficulties 
when trying to diagonalize a large nonhermitian matrix 
H' = H + W or solve the linear system (E - H')<j>{E) = 
at many values of E using general iterative techniques for 
nonhermitian matrices. For instance, the Krylov subspace 
algorithms, such as the Lanczos diagonalization procedure, 
usually converge well for the extreme eigenvalues of H', 
while numerical problems may occur for interior complex 
eigenvalues of H' in the dense spectral regions, requir- 
ing more sophisticated schemes. As such a new technique 
called PIST was recently introduced by Poirier and Car- 
rington [13]. The PIST method is based on a very efficient 
preconditioning within a QMR-Lanczos framework for it- 
erative diagonalization of large and sparse DVR Hamilto- 
nians with complex absorbing potentials. It is also appro- 
priate to mention the time-dependent approach based on 
solution of the time-dependent Schrodinger equation, 

m = e- HH 'm, (4) 

which is widely used because of its simplicity, and can 
also be viewed as a Krylov subspace method with Krylov 
vectors 

m = u i m (5) 

generated by the powers of the evolution operator U = 
e -iH DOunc j s t a t e eigenvalues = e~ tEk of U ap- 

pear at the unit circle; the resonance eigenvalues near the 
unit circle, |Afc| ~ 1 , and satisfy |Afc| < 1. That is, for a 
general initial state 0(0) all the physically important states 



tpk significantly contribute to 0(i), because they all corre- 
spond to the extreme eigenvalues of U (for which the rel- 
ative weights defined by |Afc|* are significant). This makes 
<j>(t) a convenient basis for performing the spectral analysis 
of H' . For example, one can generate a time-correlation 
function 

C(t) := T 0(t). (6) 
Because the time signal C (t) satisfies the form 

K 

C(t) :=£d fc A* fc (7) 
fe=i 

with the weights dk = Tp T Tpk'4>k4>i one can use Fourier 
spectral analysis to extract the desired eigenvalues Ek- 
Furthermore, the recently developed Filter Diagonaliza- 
tion Method (FDM) [14, 2] (see also refs. [15, 16] on 
other related superresolution methods of spectral analy- 
sis) to solve the harmonic inversion problem (7) with the 
unknowns {dk, Xk} generally leads to an enormous resolu- 
tion enhancement, thus significantly reducing the required 
propagation time, as well as the overall numerical work. 
At first glance the time-dependent framework seems nearly 
optimal as the time correlation functions can be generated 
at low cost by various iterative techniques, e.g., the split- 
operator method [17]. The difficulty, however, arises when 
both the density and the number of states are very high, 
in which case the time domain data C (t) must be very ac- 
curate at very long times in order to provide the adequate 
conditions for, e.g., FDM. Unfortunately, this requirement 
is very hard to satisfy, as it is difficult to accurately eval- 
uate the matrix exponential e~ ltH for a non-hermitian 
operator H' at very large values of t. 

Apparently, for a quantum system with Hamiltonian op- 
erator that is not explicitly time-dependent, there is noth- 
ing special about the time-dependent Schrodinger equa- 
tion, except that it provides a convenient framework for 
both thinking and devising various numerical techniques 
to solve the time-independent problem, for example, those 
based on processing the time correlation functions. As 
such we can consider alternative dynamical schemes hav- 
ing the convenient structure of Eqs. 5-7, albeit with a 
pseudo-evolution operator U, somehow related to the un- 
derlying Hamiltonian H, but whose action on a vector can 
be evaluated easily. To this end the analogy between the 
standard time evolution and the Chebyshev recursion , 

4>(t) = 2H<j)(t - 1) - 4>(t - 2) (8) 

with initial conditions 0(0) = <fr, 0(1) = H<f>, or, more 
generally, damped Chebyshev recursion [1], 

4>{t) = 2DH<t>{t - 1) - D 2 (f>{t - 2) (9) 

with initial conditions 0(0) = 0, 0(1) = DHcj) and damp- 
ing operator D (D(r) = 1 inside and D(r) < 1, outside), 
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was noticed and explored previously [18, 19, 2, 20]. In ref. 
[5] it was shown explicitly that a variant of Eq. 9, 

4>(t) = 2DH(f>(t-l)-D<f>(t-2) (10) 

with initial conditions <fi(0) = <j>, 0(1) = 0, can be writ- 
ten in the familiar form (5) with some effective evolution 
operator U. The corresponding pseudo-time Schrddinger 
equation allows one to reap all the benefits of the time- 
dependent methods without having to deal with the stan- 
dard time-dependent Schrodinger equation (4) involving 
nonhcrmitian Hamiltonian; it only requires the evalua- 
tion of a single iJ-matrix-vector product per time step 
and avoids the use of complex arithmetic, even when the 
absorbing potential is implemented. 

As is well known, a (physical time) autocorrelation func- 
tion at time 2t can be computed by solving the time- 
dependent Schrodinger equation up to time t, since one 
can use 

C(t) :=<f> T U*<t>=(U t <f>) T (U t <f>), 

assuming only the complex symmetry of the evolution op- 
erator U, which is always easy to achieve, both with ab- 
sorbing potential or without. For the Chebyshev autocor- 
relation function 

y(t) := <j> T <f>(t), (11) 

where the vectors cj>(t) are generated by Eq. 8, a factor of 
two saving is also well known (see, e.g., the discussion in 
ref. [2]): 

y(2t + p) = 2cf ) (t) T <t > {t+ P )-y(p) 

with p = 0,1. When the damped scheme (9) is imple- 
mented, this recipe does not provide the correct autocor- 
relation function as defined by Eq. 1 1 . Nevertheless, it was 
tried by Li and Guo [4]. The resulting doubled sequence, 
when processed by FDM, gave approximately correct reso- 
nance energies and widths. In ref. [5] starting with the re- 
cursion formula (10) we derived an exact doubling scheme 
for the corresponding autocorrelation function (11): 

y(2t + p) = <j){tf >(t + p) - <f>(t + l) T D"V(t + 1 + P) 

with p = 0, 1. 

We note that the damping operator D in Eq. 9, as well 
as in Eq. 10, effectively leads to an energy-dependent com- 
plex absorbing potential We{t), which may have a small 
or large real part relative to its imaginary part, depend- 
ing on the particular recursion formula and the energy 
E. Usually, having too large Re We{t) is not desirable. 
Unfortunately, within the two frameworks one has little 
control over this circumstance. In the present paper, we 
re-derive and extend the above results by using a more flex- 
ible form for the absorbing potential, which can be adapted 
easily. The resulting pseudo-time Schrodinger equation is 
shown to be very convenient for calculating various dy- 
namical properties, such as resonance parameters or ma- 
trix elements of the Green's function. The latter is carried 



out by generating the pseudo-time cross-correlation func- 
tions followed by their inversion. A numerical example 
demonstrates the validity of the theory and compares dif- 
ferent damping schemes, suggesting an optimal scheme, 
that leads to accurate results with minimal computational 
effort. 

QUADRATIC EIGENVALUE PROBLEM 
FROM NONLINEAR SPECTRAL MAPPING 

From now on, we assume that the Hilbert space is dis- 
cretized so that the states are vectors ip € C K and H, W 
are real symmetric K x K matrices, W diagonal, as, e.g., 
in the case of a discrete variable representation [27]. 

Let (A)$ := tp* Aip/ip*ip define the expectation value 
of the operator A. The original Hamiltonian matrix has 
its spectrum between H — AH and H + AH. The quan- 
tity AH is called the spectral range and is defined by the 
particular basis set used to represent the Hamiltonian. 
However, to simplify the following equations we assume 
without loss of generality that the Hamiltonian matrix is 
already shifted and scaled so that 

\(H-D )^\<1 (12) 

for any state ip, where the diagonal real symmetric oper- 
ator Da will be specified later. (In places where the scale 
factor AH is relevant, it will be inserted explicitly.) 

Such a scaling is implemented routinely in the frame- 
work of the Chebyshev polynomial expansion. We note 
that this typically moves the ground state energy to the 
lower edge of the spectrum, E = — 1, while the energies 
of interest, including the dissociation energies, appear at 
the bottom of the spectrum, E <~ —1, as the total spec- 
tral range is usually an order of magnitude larger than the 
physically relevant spectral range. 

We consider the choice 

E-W = D n + ^uD 1 + \- 1 D 2 , (13) 

for some complex parameter u and real matrices Do, D\ 
and D 2 , diagonal in the coordinate representation and sat- 
isfying 

D (r) = 0,Di(r) = l,D 2 (r) = l for r e ft. (14) 

This is a generalization of the previous results where the 
special cases using D = 0, D 2 = Df 1 (as in ref. [1]) and 
Dq = 0, D 2 = 1 (as in ref. [5]) were encountered. 

To match the original problem in ft, where the absorbing 
potential W(r) vanishes, u = ue and E = E(u) must be 
related by 



u E = E + Wl-E 2 , or E(u) = -u+-u- 1 (15) 
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as shown in Fig. 1. The absorbing potential then becomes 
a function of energy E (or, equivalently, of u = ue)- 

W B = ~u(l-D 1 ) + ~u- 1 (l-D 2 )-D a (16) 

Insertion of (13) into (2) gives a nonlinear eigenvalue prob- 
lem for it, 

= (d + X -uD x + \u- X D^ ^. (17) 
We may think of this equation as an eigenvalue problem 
Hip = E{u)ip (18) 
involving an operator-valued it-dependent energy 

(19) 



E{u) =D + ^uD 1 + \u~ 1 D 2 



that, by (14), reduces in fl to the constant E = ^u+^u 1 




RECASTING TO AN ORDINARY 
EIGENVALUE PROBLEM 

Apparently, Eq. 17 is equivalent to a linear eigenvalue 
problem in a space of doubled dimension: 



Utp = wtp 



(20) 



with 2K dimensional state vectors tp = (X) and square 



2K x 2K matrix 
U := 



/ 



D^D 2 2D^ 1 (H-D )J' 1211 

where / denotes the K x K unit matrix. This can be seen 
since (20) yields 

ip' = uip, D^D 2 ijj - 2uDi\H - D )ip + u 2 i) = 0. 

Now multiplying by Di/2u, we find (17). 

We have rewritten the original nonhcrmitian eigenvalue 
problem (2) as another nonhcrmitian eigenvalue problem 
(20), but with the matrix U of doubled dimension. This 
would not necessarily be an advantage, if the eigenvalues 
of U did not possess the very important property that the 
spectral domain of U is the unit disk. 

To see this, we consider an eigenpair (u,ip) of (17) with 
tp*ip 7^ 0. Multiplying (17) by 2uip* gives the quadratic 
equation 

u 2 (Di)^ - 2u(H - Dq)^ + (D 2 ) t = 0. 
with solutions 

(H - DoU ± i^J {Di)^{D 2 )^ — (H — D ) 



(Di) 



(22) 



To guarantee correct behavior in the following results 
we need to apply the componentwise restrictions on the 
diagonal damping operators: 



< Dl 1 < D 2 < -Di > 1, 



(23) 



FIG. 1. The spectral mapping Eq. 15 maps a bound state E\, 
which is real, to the upper half of the unit circle, and the resonance 
states, E2 and £3, with Im < 0, inside the upper half of the unit 
disk. «2 ~ 1 as it corresponds to a narrow resonance. (In practice 
the physically relevant states appear in a small subset of the whole 
spectral domain.) 

which together with the condition (12) imply that the 
square root is real. Thus, the solutions of (17) come in 
complex conjugate pairs, (it, ip) and (u, ip). The physically 
relevant eigenenergies with Im E < come from u with 
Im it > 0. (22) implies that 



(H 



D )l + {D 1 )^{D 2 U-{H-D ) 



) 2 

It), 



(D 



Pi)* 



< 1 



(24) 



by (23). Thus, u is a complex number lying in the upper 
half of the unit disk (see Fig. 1). Moreover, |u| = 1 if and 
only if (W)-ip — 0, i.e., if and only if ip has support in il, 
which is the case for the bound states. The states with 
|u| ~ 1 correspond to the narrow resonances. 

The eigenpairs {uk,ipk) of (17) can be used to evaluate 
the physically interesting quantities (e.g., the complex res- 
onance energies Ek, scattering amplitudes, etc.). However, 
because of the nonlinearity they have somewhat different 
properties from those of the regular nonhermitian eigen- 
value problem, which we now proceed to derive. 

COMPLETENESS 

As was already established, the eigenpairs (itfe, ipk) of U 
satisfy 



4>k 



(25) 
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where (u/ s ,V'fe) is an eigenpair of (17). Since conversely, 
any such eigenpair determines an eigenpair of U, the 
nonlinear eigenvalue problem (17) has at most 2K dis- 
tinct eigenvalues. If there are 2K distinct eigenvalues 
Mi, ... , U2K, the matrix U is diagonalizable, and there is a 
basis x/ji, . . . , ^2K of eigenvectors of the form (25). There- 
fore, we can obtain the completeness relations, i.e., for any 
vector (f> we may write 



2K 



2K 



k=\ 



Oku k ipk 



(26) 



with uniquely determined coefficients 9 k . 

ORTHOGONALITY 

Using (17), the symmetry of D , D\, D 2 and H, and the 
fact that Di commutes with D 2 , we may compute xpj Hip k 
in two different ways: 



^]H^ k = Vj (D + -M feJ Di + 



Do 



2«' 



D' 



ipk 
ipk- 



For j k, assuming that the eigenvalues are not degen- 
erate, we may take the difference, multiply by the factor 
2ujUk/(uk — Uj), and find that the eigenfunctions satisfy 
the orthogonality relations 



ipj (D 2 - UjU k Di) ip k = 0. 



(27) 



Provided that the left hand side of (27) does not vanish 
for j = k, we may normalize the eigenvectors so that 



tpj {D 2 - UjUkDi) ip k = S jk . 



(28) 



Due to (25), the orthogonality relations can be rewritten 
in the double-dimension form: 



^p]D%p k = 5 jk with D := 
Using (29) and (26) we find: 



D 2 








(29) 



9 k = ^OjSjk = Y,fii$jf*i>k = $ T D^ k (30) 



^2u k ip k tp k = 0, 



fc=i 

2K 

£ 

k=i 



ulipkipl 



Di 



THE GREEN'S FUNCTION 

For ue = E + i\J\ — E 2 and an eigenpair (u k ,tp k ) of 
(17) we can write 



(H-E + W E )^k 



UE - Uk 
2UkUE 



(D 2 - UkU E Di)ip k - 



Now multiplying from the right by 



2u k u E 
u E -u k 



v> fe T , 



summing 



over k and using the relations (32) we obtain 

2K 

{E-H-We)^^^-^ 



k=l 



U E - U k 



2K 



^2(D 2 - u k u E Di)tp k tp k = I. 



k=i 



This leads to a spectral representation of the damped 
Green's function G W {E) := (E - H - We)' 1 



(33) 



This expression uses the eigenvectors, which may be diffi- 
cult to generate in large-scale computations. However, we 
shall see that (33) leads to iterative schemes for computing 
matrix elements of Gw{E). To derive these, we introduce 
the double-dimension Green's function 



k—1 



G W (E) * 

* * 



(34) 



whose top left corner suffices to represent the matrix el- 
ement of the Green's function. Indeed, for suitable initial 
states 



and final states <p a , we have 



4>lG w {E)<p = 4>lG w (E)fo 



(35) 



RESOLUTION OF IDENTITY 

The orthogonality relations (29) and completeness imply 
the resolution of identity, 



2K 



2K 



J2^Id = dJ2^1 



k=l 

which in the explicit form reads 

2K 



k=l 



D 



2 ! 



(31) 



(32) 



with 











Gw{E) is representablc in terms of U as 
G w (E) = {l-U/u E )- 1 Ub-\ 



(36) 



(37) 



k=i 



which can be verified by applying (34) and (37) to Dip k . 
(Note that the term D^ 1 in (37) is unnecessary as one 
usually considers initial states with support in f2, where 
D^ 1 = 1, but we still prefer to keep it for completeness.) 
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THE PSEUDO-TIME SCHRODINGER 
EQUATION 

By replacing (1 — U/ue)^ 1 in (37) with the geometric 
series we can express Gw{E) as a power series in U : 



G w (E) = J2u t E 1 U t D- 
t=i 



(38) 



This form using a discrete Fourier transform of the pseudo- 
evolution operator 17* is reminiscent of the integral Fourier 
transform of the true evolution operator 



G(E) 



•f 

Jo 



-iHtAEt 



dt. 



Howc\ 



cr, for a given absorbing potential and basis set, 
can only be represented approximately while Eq. 
38 is exact as there is no approximation involved in eval- 
uating U . Because U is bounded inside the unit disk (cf. 
Eq. 24), the expansion (38) also leads to a numerically 
stable procedure to compute the Green's function matrix 
elements 



]G w {E)4>p 



1 y a p{t), 



(39) 



where we have introduced the pseudo-time correlation 
function 

y a[i {t) := <f>ZMt) = &pit) (40) 
= ^U'D- 1 ^ (t = 0,l,...). 

Eqs. 39 and 40 constitute an important result as they 
provide a convenient and efficient numerical framework 
based on solving the pseudo-time Schrodinger equation 

^i^^ii'b- 1 ^ (t = o,i...). (4i) 

That is, the quantum dynamics problem is finally reduced 
to the Fourier spectral analysis of pseudo-time correlation 
functions. 

By noticing that the state vectors 4>p(t) have the form 

(41) can be rewritten as a variant of a damped Chebyshev 
recursion formula 

<f>p(t) = D^ 1 [2(77 - D )Mt ~ !) - D Mt - 2)] , (42) 

£ = 2,3,... 

with initial conditions 4>fj(0) = D^&p and 0/3(1) = 0. 

The numerical cost of implementing a single step in (42) 
is essentially equal to that of multiplication of a real sym- 
metric K x K matrix H by a real vector, as the other 
matrices are diagonal. 

TIME DOUBLING 

We note that DU* is a complex symmetric matrix, 
(7)t7*) T = (C/*) T £> = DU\ 



which follows from the spectral representation 

IK 



(43) 



fe=i 



or can be verified directly by the matrix multiplication. 
Thus, for any t and s we can write 

(44) 



where <j) a (s) and 4>p(t) are the solutions of Eq. 41 with 
initial conditions Q (O) = 7) _1 Q and (f)p(0) = D^ 1 ^. 
This means that y a p(t) can be computed up to time 2T 
using 

y al3 (2t+p) = 4> a {t) T DMt+p), p = 0,l, (45) 

concurrently with the computation of <f> a {t) and 4>p(t) for 
t = 0, . . . , T. In particular, the calculation of an auto- 
correlation function y aa {t) up to time 2T requires ~ T 
multiplications of the real and sparse K x K i7-matrix on 
a vector with just a few vectors stored at a time. 

RESONANCE CALCULATION BY 
HARMONIC INVERSION OF PSEUDO-TIME 
CORRELATION FUNCTIONS 

The spectral mapping Ek — ► Uk (15) moves all the phys- 
ically relevant eigenvalues to the vicinity of the unit circle, 
i.e., they all become extreme (or nearly extreme for the 
resonances) eigenvalues of U. This, in turn, creates very 
favorable conditions for computing these eigenvalues iter- 
atively using any suitable Krylov subspace method with 
the Krylov vectors generated by the powers of U, because 
the extreme eigenstates contribute most to the latter. In 
the present framework, such a strategy of computing the 
bound and resonance state energies, in its most simple and 
convenient form, boils down to the harmonic inversion of 
pseudo-time correlation functions (40), which due to (43) 
satisfy 



2K 



y a f){t) =^2d a0k u t k 



(46) 



k=i 



with 



dafik — b a k 



b k = (J>li>k) (i>k<i>p) ■ 



(Note, that the resolution of identity (31) implies 

d Efc V'fcV'fe 4> a = DJ2k Kk^k-) 
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Thus the nonlinear eigenvalue problem (17) is reduced 
to the signal processing problem of finding the spectral 
parameters (u k ,d a ^ k ) (k = l,...,2K) satisfying (46) 
for the sequence(s) y a p(t) computed by (42) and (40). 
The simplest procedure then corresponds to considering 
a single doubled autocorrelation function y aa {t) with t = 
0, . . . , 2T. In exact arithmetic the harmonic inversion of 
such a sequence will give the exact results if T > 2K, 
thus, using only T ~ 2K of real matrix-vector products. 
However, this is impractical as it would formally require 
to solve aTxT eigenvalue problem. To reduce the com- 
putational burden and to maintain numerical stability the 
eigenvalues are extracted very efficiently in a small Fourier 
subspace by the FDM [14, 2]. If we introduce the unsealed 
energy e — EAH + H, the required length 2T of the dou- 
bled sequence needed to converge an eigenenergy E k (cf. 
Eq. 15) by the FDM will be defined by the locally aver- 
aged density of states p{e) for E k ~ E and the spectral 
range AH of the Hamiltonian matrix according to the ap- 
proximate relationship [2, 31] 

T > 27rA#p(£)Vl -E 2 , (47) 

where the factor y/l — E 2 arises from the u — > E mapping 
(15). 

Because the eigenpairs come in complex conjugate pairs 
(u k ,ip k ) and (uk,4>k), once the initial vectors (\> a and (f)p 
are real (46) becomes 

K 

y a p{t) = Re d <*f3ku{, (48) 
fc=i 

where only the physical eigenvalues with Im u k > are 
included in the sum. Note also that for the pure bound 
state problem, when \u k \ = 1 and the eigenfunctions tp k 
(and, therefore, the coefficients d a p k ) are real, the se- 
quence y a p{t) has the time reversal symmetry 

y a p{-t) = y a p{t) (49) 

which further doubles the total time by extending the sig- 
nal to the negative times. 

INVERSION OF THE TIME 
CROSS-CORRELATION FUNCTIONS 

As was previously shown in rcf. [3] one can, in prin- 
ciple, compute matrix elements (fi^Gw(E)(j>0 for any (f> a , 
4>p by (i) propagating a single initial state 4>o using (42), 
(ii) computing the cross-correlation functions yoo(t), Voa(t) 
and 2/o/3(0 using (40) and (iii) solving the corresponding 
harmonic inversion problems (46) for the unknown param- 
eters {life, b akl bfj k } to evaluate 

<f>lG w (E)cP = V ^h^- bak b fik . (50) 



Moreover, a more stable evaluation of the matrix element 
( t > H l Gw{E)(f)p may be achieved using the Regularized Re- 
solvent Transform (RRT) which was described in detail in 
ref. [15]. The advantage of RRT is that the calculation of 
spectral parameters {u k , b ak , bp k } is avoided, the spectra 
are computed directly by matrix inversion or by solving 
linear systems. 

This results in an enormous numerical saving as, tradi- 
tionally, the time-dependent approaches require multiple 
initial state propagations in order to compute, for example, 
the full S-matrix or cumulative reaction probability. How- 
ever, the described approach, although formally exact (in 
exact arithmetic), is applicable only when the dynamics is 
governed solely by narrow resonances: the parameters of 
very broad resonances (or poles of the Green's function) 
are generally grossly inaccurate, leading to very unstable 
spectral estimation by Eq. 50. A significant numerical 
saving in cases with broad resonances is, however, achiev- 
able if one adapts another strategy in which one (i) prop- 
agates a set of initial states {4>fj}, {fi = 1, ■ • • ,L) using 
Eq. 42, (ii) computes the cross- correlation matrix using 
the doubling formula (45) and (iii) solves the harmonic 
inversion problem (46) for the L x L x 2T data matrix 
(see refs. [14, 3, 15]). As argued in rcf. [3], this approach 
has a potential of reducing the total propagation time T 
required for the accurate harmonic inversion by a factor 
of L, i.e., preserving the total number of matrix-vector 
multiplications, N tota \ = L x T. Unfortunately, to use the 
doubling trick (45) in this case, the storage requirement to 
generate the cross-correlation matrix is increased because 
of the need to simultaneously propagate L states rather 
than one. 

PRACTICAL CONSIDERATIONS: CHOOSING 
THE ABSORBING POTENTIAL 

The physical observables can generally be computed 
from the Green's function matrix elements at real ener- 
gies. For this reason and for the sake of simplicity in the 
following analysis we will assume the energy E to be real. 
However, essentially the same conclusions hold for com- 
plex energies near the real axis, i.e., including the narrow 
resonances. 

Since u E x = E — i\f\ — E 2 , our construction (16) leads 
to the complex valued absorbing potential 

W E = E I 1 I - A) -«Vl -E 2 

(.51) 

The restrictions (23) on the choice of D\ and D2, im- 
ply that Im We < 0. This is essential since it leads to 
the correct limit for the Green's function Gw{E) on the 
real line. The real part Re We is generally nonzero; al- 
though its role is not crucial in scattering calculations, 
significant values of Re We may or may not be desirable. 
In particular, a positive Re We reduces the total density 
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of states p(E), while affecting only slightly the resonance 
states, and thus, it may accelerate the convergence (cf. Eq. 
47). At the same time Re We may get excessively large 
compared to Im We, for example, near the spectral edge, 
E <~ — 1, if Do = 0. It is not absolutely clear without nu- 
merical tests how this artifact will affect the results, while 
it can be controlled if the matrix Do, which is in principle 
unrestricted, is used. Namely, relatively to the total spec- 
tral range AH (which becomes unity after scaling of H), 
the physical energy range is generally only a small part of 
the unit disk. Therefore, within this range we can often 
assume We ~ We for some reference energy En. Now by 
setting 

D 1 + D 2 s 



D Q = En 1 - 



(52) 



we can significantly reduce (or increase) the real compo- 
nent if needed, 



Re W E = (E - E) [ Dl+D2 _ - 



(53) 



without affecting Im We- In order to avoid negative 
Re We in the physically relevant energy region, it suffices 
to take Eo higher than the maximum energy of interest. 

Clearly, the behavior of the absorbing potential We 
where it starts to turn on has the main effect for the scat- 
tering calculations. In fact, as follows from our numerical 
tests, the best performance is achieved when both D\ k> 1 
and Di « 1 over almost the entire grid used to represent 
the state vectors. Therefore, it is convenient to use the 
real operators 71 and 72 defined by 



D 2 



(54) 



which vanish in f2 and slowly turn on in the absorbing 
region, and assume that both 71 and 72 are small. By 
expanding into the Taylor series up to the second order, 
we can rewrite (51) as 



W E « (Eq-E) 



- Wl-E 2 



7i - 72 
2 

71 + 72 



7? 



7l 



(55) 



+ 



7i 2 



7l 



We expect the results to be generally insensitive to a 
particular form of We as long as 71 and 72 are sufficiently 
smooth and have sufficiently large spatial extension. So 
here we only consider the two special choices [1, 5] 

M&T scheme: 71 = 72 = 7 as in ref [1] with j(R) > 0: 

W E = (Eo - £)(cosh( 7 ) - 1) - iyjl - E 2 sinh( 7 ).(56) 

As follows from the expansion (55), Re We is already small 
as its leading term becomes quadratic in 7, while Im We 
is linear: 

We « \(E - £)7 2 - iy/l - E 2 -f. (57) 



Note again that in the computations involving molecu- 
lar vibrational spectra the energies of interest are usually 
close to the bottom of the spectral range of the scaled 
Hamiltonian matrix, E <~ — 1. Defining the shifted energy 
s := E + 1 with zero at the bottom of the spectrum and 
assuming £ < 1, we can approximately determine how the 
absorbing potential depends on energy: 



W E ~ ^(eo - e)7 2 - iV2ej- 



(58) 



A convenient choice for 7 in Eq. 56 corresponds to ■y(R) 
being a function of the reaction coordinate R, vanishing 
in the interaction region, R < Rn, and smoothly growing 
in the absorbing region, Rq < R < -R max , where i? max 
defines the fartherst grid point: 



70R) 



A 



R — Rn 
AH \R max Rn 



(59) 



Here A is an adjusting parameter. As follows from Eq. 
58, the factor v AH in Eq. 59 minimizes the sensitivity of 
Im We to the actual spectral range AH of the Hamiltonian 
matrix. Note that with this construction Re We will be 
sensitive to AH, but generally small. 

N&M scheme: 72 = and 71 = 27, as in ref [5]: 
Eq-E- Wl - E 2 



W E = 



(e 2 -> - 1) (60) 



(E -E- iVl-E 2 )~f, 



(61) 



where we retained only the leading linear term. Further 
assuming e to be small, we obtain 



W E ~ (eo - e)l - iV^ej- 



(62) 



Thus the difference between M&T and N&M schemes is 
the different behavior of Re We ■ This difference may be 
eliminated almost completely by manipulating with Do- In 
the next section, we examine the effect of this difference 
on the quality of the resonance calculations. 

In order to obtain the best resonance energy estimates 
for given basis set, the common practice is to search for 
stationary eigenvalues Ek under variations of the absorb- 
ing potential [7, 10, 23]. This is usually done by varying 
its amplitude, which here corresponds to the free parame- 
ter A. (One could possibly change instead Eo, or treat Eo 
as a function of A.) The eigenvalue trajectories Ek(X) are 
then analyzed for stationary points, such as the point of 
maximum curvature, which in an ideal case may be a cusp 
(see Fig. 2). One way to practically find such a point is 
to minimize the derivative dEk/dX. The generation of the 
eigenvalue trajectories usually significantly increases the 
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computational time, while a reasonable accuracy for most 
states can be achieved from a single or a few runs using 
an a priori established optimal value of A. Therefore, such 
a search is justified only if a very high accuracy for all 
resonances is needed. The cases where eigenvalue trajec- 
tory analysis may be necessary include extremely narrow 
(shape) resonances. 

NUMERICAL EXAMPLE: RESONANCES OF 
THE HCO MOLECULE. 

In this section we apply the present approach for an ac- 
curate computation of the nonrotating HCO resonances. 
This system has been used in the past by various groups 
as a benchmark problem to test new approaches (see the 
comprehensive review [13] and references therein). In par- 
ticular, the first resonance calculation using FDM was ap- 
plied to HCO [24]. We use the potential energy surface 
of Keller et al [25], where a resonance calculation was also 
performed and the resonances were assigned. Recently, 
Poirier & Carrington [13] applied the above mentioned 
PIST method to reproduce these results but with a sig- 
nificantly higher accuracy and numerical efficiency. We 
use the latter results as a reference. Although, PIST is a 
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FIG. 2. Typical eigenvalue trajectories for two resonance states 
obtained using three different damping schemes (see text). 



very efficient alternative to the present technique, in our 
study we do not make an adequate comparison of the two 
methods as it would require the use of, at least, the same 
basis set. 

The present choice of the basis is the most primitive 
direct-product Jacoby coordinate DVR grid, similar to 
that used in ref. [24]. Several basis sets with various grid 
sizes and cutoff parameters were tested. The results re- 
ported below correspond to a particular single set of pa- 
rameters, which, as was verified empirically, provided ex- 
tremely accurate results for the resonances in the energy 
region below 9000 cm -1 . In order to avoid an extensive 
search for optimal basis set parameters, our strategy is 
to use a larger and denser grid than may be necessary to 
garantee high accuracy. This strategy is justified by a very 
favorable numerical scaling of the present technique with 
the basis size. 

For i?H-co (the dissociation coordinate) we used 160 
sinc-DVR [21] points in the interval [2,8] ao, contracted 
to N R — 40 points by means of the HEG [22] method 
using the eigenfunctions of the ID Hamiltonian defined 
by V(R) = min r , e V(i? e ,r,0), where V(R e ,r,6) is the 3D 
potential of HCO in the Jacoby coocrdinates. For the 
rco coordinate we used 64 sinc-DVR points in the interval 
[1.8,3.5] ao, contracted to N r — 16 points. For the angu- 
lar variable we used Ng — 46 Gauss-Legendre-quadrature 
DVR points. In order to reduce the Hamiltonian spectral 
range AH we replaced the high values of the potential 
by the cutoff value V cut = 25000 cm -1 . This resulted in 
AH = 44654 cm -1 . We note, that for the present ap- 
proach the quality of the grid basis is measured by AH, 
which directly affects the convergence (cf. Eq. 47). Since 
our basis is very primitive the spectral range is relatively 
high. One may argue though that a primitive basis has 
the advantage of minimizing the number of adjusting pa- 
rameters and leads to a fast matrix-vector multiplication 
for given grid size. 

Here we report results corresponding to the three dif- 
ferent schemes: (I) N&M scheme with Do defined by 
Eq = 9000 cm -1 (the maximum energy of interest), (II) 
N&M scheme with D = 0, and (III) M&T scheme with 
Dq = 0. (We note that the use of nonzero Do in the M&T 
scheme has almost no effect as Re We is already very small 
in this case.) 

For each scheme the damping potential has been taken 
in the form (59) with Rq = 5 ao- For this study we varied 
the strength parameter A in the range [0.0002, 0.6] with a 
logarithmic distribution of 22 A values. 

For each A an autocorrelation function y(t) (t = 
0, 2T — 20000) with a random initial vector was gener- 
ated using T = 10000 matrix-vector multiplications. Cal- 
culation of a single autocorrelation function takes 239 sec 
on an Athlon 1.8 GHz pc using the gnu g77 Fortran com- 
piler. Because the latter is very inefficient, the reported 
timing only gives a rough estimate of the algorithm effi- 
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ciency. The autocorrelation functions were then inverted 
by our quadruple-precision FDM code [26]. For a large 
enough value of A, the FDM results converge for most res- 
onance states using T w 7000, however, to maintain an 
extremely high accuracy, especially for the lowest sharp 
resonance |013) with width T = 3 • 10~ 8 cm -1 and for 
small A values, we needed T 10000. (The effective den- 
sity of states, that includes both the true resonance states 
and states representing the continuum, becomes high for 
small A.) Of course, a smaller grid and lower V^ut would 
result in a smaller AH reducing the needed number of it- 
erations accordingly. Also note that because only small 
generalized eigenvalue problems are encountered in FDM, 
the cpu-time increase due to the use of the quadruple pre- 
cision is not essential, while it improves the accuracy and 
accelerates the convergence of the extracted eigenvalues 
significantly, compared to the double precision code. For 
the present "toy problem" the harmonic inversion part is 
as time-consuming (with quadruple precision) as the cor- 
relation function generation, however, the former becomes 
relatively negligible for bigger systems. 

Let -E r os := Re Ek define the positions and Tk '■— 
— 21m Ek, the widths of the resonances. The eigenvalue 
trajectories Sfc(A) with Tk < 120 cm -1 are plotted in the 
complex plane and the most stationary point (see above) 
is assumed as the best resonance energy estimate. In Fig. 
2 we show the eigenvalue trajectories generated using the 
three different damping schemes decribed above for two 
states (a) with E Tes = 3999.2, T = 32.8 cm" 1 and (b) 
E rcs = 1251.148, r = 1.8 ■ 10" 6 cm" 1 . Schemes I and III 
result in quite different trajectories for most A values but 
have easily identified stationary points (cusps) for about 
the same value A. These stationary points are very close to 
each other. At the same time, scheme II gives more com- 
plex trajectories that circle around the resonance position, 
but have no well identified stationary point. Similar pat- 
terns are observed for most other states. Note also, that for 
some resonances the cusps may not exist, however having 
the two types of trajectories helps one to better locate the 
complex resonance energy by taking the point where the 
trajectories of different types approach each other. Fur- 
thermore, scheme III generally results in sharpest cusps, 
but because it corresponds to a smaller Re We, it also 
requires larger values of T for an accurate harmonic inver- 
sion, as smaller Re We corresponds to higher density of 
states. For the same reason scheme II needs the smallest 
T. 

Surprisingly, for most resonances, the stationary point 
in the eigenvalue trajectories are approached for the same 
value A opt = 0.08, while only for a few extremely sharp 
resonances near the dissociation threshold E^is — 1086 
cm -1 , the values Ek{\ = 0.08) give poor estimates of the 
true complex resonance energies. For these states the op- 
timal A is a monotonically growing function of E reB — i?dis i 
the optimal A for the sharpest state, |013), is A op t ~ 0.0016 



and for the next state, 1 005) , it is A op t ~ 0.0064. By set- 
ting A = 0.08 with the present choice of j(R), one can use 
a single autocorrelation function to obtain good resonance 
estimates for most states, but a few. 

There is no simple way to determine the error bounds 
other than by comparing the results using different pa- 
rameters. The accuracy in the positions of the resonances 
varies depending mostly on Tk- It also depends in a less 
systematic way on the quantum numbers. For narrow res- 
onances (Tk ~ 1 cm -1 or less) the accuracy in the position 
is of the order of 0.01 cm- 1 or better, while for the other 
states with Tk < 100 cm.- 1 it is generally of the order of 
0.1 cm -1 . We believe, that generally the error in the re- 
ported widths is less than 1%, although this may not be 
true for some states. 

Fig. 3 and Table I summarize our results obtained 
with scheme I and compare them with the calculation of 
Poirier & Carrington [13]. The assignments of the reso- 
nance states have been made by Keller et al [25]. The 
latter results for the resonance parameters are not dis- 
played here as they are less accurate. Our results generally 
agree very well with the results of Poirier & Carrington, 
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FIG. 3. The widths T versus resonance positions E TCB shown for 
two different scales. The present results (circles) are obtained using 
scheme I (see text). The crosses are the results of Poirier & Carring- 
ton [13]. 



PSEUDO-TIME SCHRODINGER EQUATION 



11 



although some states are missing in the latter, such as 
the broad resonance shown in Fig. 2: missing eigenvalues 
are a common drawback of Lanczos based diagonalization 
approaches. Note also that for the first two sharp reso- 
nances, Poirier & Carrington did not report reliable width 
estimates. 

Even though under the present circumstances a fair 
comparison between the present approach and PIST [13] 
is not possible (e.g., very different basis sets were used in 
the two cases), the following remarks seem appropriate. 
At least for the HCO case, PIST appeared very efficient 
(required relatively few matrix-vector multiplications to 
accurately compute the resonances) and relatively small 
cpu-time, although the fact that it has more adjusting pa- 
rameters than the present technique may be viewed as its 
disadvantage. To conclude, the present approach, at least, 
for the resonance calculations, appears very reliable, accu- 
rate and efficient. It requires a minimal number of adjust- 
ing parameters and scales favorably with the size of the 
system. The method is quite flexible in the choice of the 
damping scheme, but according to our tests an optimal 
choice corresponds to scheme I (see above) with A w 0.08. 

We make the corresponding Fortran codes available 
upon request (e-mail: mandelsh@uci.edu). 
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TABLE 1 

Resonance energies E reB and widths T (cm -1 ) for the nonro- 
tating HCO molecule. The quantum numbers describe, respec- 
tively, the CH stretch, the CO stretch, and the bend. (The 
last digit in the "present results" gives a rough estimate of 
the computational error, although an accurate error estimate 
is generally unknown.) 
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3409.6 


17.17 


3409.66 


17.18 


052) 


7233.1 


18.2 


7233.23 


18.26 


211) 


3512.31 


22.43 


3512.29 


22.39 




7281 


74 


7281.1 


74.1 


105) 


3566.1 


57.8 


3565.61 


57.9 


|044) 


7493.83 


8.73 


7493.84 


8.68 


032) 


3704.872 


0.238 


3704.866 


0.240 




7568 


50.3 


7563.4 


52.2 


024) 


3921.77 


4.677 


3921.77 


4.687 




7585.8 


74.8 


7587.2 


68.6 


130) 


3999.9 


18.9 


3999.84 


19.80 


150) 


7602.67 


7.31 


7602.63 


7.29 








3999.2 


32.8 


036) 


7656.7 


29.56 


7656.64 


29.68 


|016) 


4036.49 


10.03 


4036.63 


9.94 


142) 


7813.3 


47 


7813.1 


48.1 


220) 


4084.6 


19.5 


4084.71 


20.25 


240) 


7884.9 


66.6 


7884.9 


66.2 


122) 


4214.7 


23.4 


4214.49 


23.12 


061) 


7952.3 


5.53 


7952.29 


5.52 


114) 


4345.8 


23.68 


4345.79 


23.80 






8030.7 


74.3 








4392.1 


113.4 


224) 


8171.1 


105 


8170.6 


105.1 


|041) 


4436.7 


9.9 


4436.61 


10.03 


151) 


8232 


30.3 


8233.06 


29.48 




4463.8 


12.45 


4463.70 


12.44 


053) 


8237.3 


94 


8236.6 


96.1 


|212) 


4514.2 


47.65 


4514.1 


47.8 




8287.1 


68.2 


8287.2 


66.9 


033) 


4726.773 


0.728 


4725.768 


0.727 


|045) 


8445.08 


37.92 


8444.9 


37.9 


017) 


4885.6 


10.4 


4885.78 


9.98 




8488 


41 


8486.7 


37.5 








4887.8 


65.6 




8556.2 


53 


8557.5 


52.6 








4938.4 


112.9 


151) 


8558 


38 


8557.7 


39.0 


|131) 


5057.9 


17 


5057.81 


17.12 


070) 


8616.5 


0.728 


8616.421 


0.717 


213) 


5143.5 


73.2 


5143.6 


73.5 


241) 


8684 


77 


8687.8 


76.3 


050) 


5156.2 


0.054 


5155.939 


0.0544 


143) 


8802.3 


90 


8801.8 


88.6 


123) 


5275.7 


32.1 


5275.8 


32.3 




8921.4 


90.9 


8921.4 


90.9 


221) 


5400.2 


46.8 


5400.5 
5460.9 


46.6 
114.4 


|062) 


8973.7 


7.7 


8973.27 


8.28 



